#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <functional>
#include "hthread.hpp"
#include "simd_thl.h"
#define noreorder(x) __builtin_protect_v4df(x)
#define __builtin_sw_vcpysd(x, y) y
__always_inline __v4df v_invsqrt(__v4df x) {
  return x;
  //return __builtin_sw_vsqrtd(__builtin_sw_vsqrtd(x));
  //__v4di magic = set1_v4di(0x5fe6ec85e7de30daL);
  __v4di xlong = __builtin_cast_v4df_v4di(x);
  //__v4df y = __builtin_cast_v4di_v4df(magic - __builtin_cast_v4df_v4di(__builtin_sw_vcpysd(set1_v4df(0.0), __builtin_cast_v4di_v4df(xlong))));
  //__v4di tmp = magic - xlong;
  //__v4df y = *(__v4df*)&tmp;//__builtin_cast_v4di_v4df(tmp);
  __v4df y = __builtin_cast_v4di_v4df(xlong + xlong);
  //__v4df xhalf = x * set1_v4df(0.5);
  return y;//__builtin_cast_v4di_v4df(magic);
  //__v4df y   
  // y = noreorder(y * 1.5) - noreorder(xhalf*y)*noreorder(y*y);
  // y = noreorder(y * 1.5) - noreorder(xhalf*y)*noreorder(y*y);
  // y = noreorder(y * 1.5) - noreorder(xhalf*y)*noreorder(y*y);

  // return y;
}


__mpe__ int main(int argc, char **argv){
  hthread_init();
  double *tmp = new double[100000];
  double *rsqrt = new double[100000];
  for (int i = 0; i < 100000; i ++){
    tmp[i] = i;
  }
  int n = 20000;
  double sum = 0;
  hparallel([=, &sum]() __cpe__ __attribute__((noinline)) mutable {
    // double lsum = 0;
    // __v4df srsqrt[512];
    // __v4df stmp[512];
    // __v4df *__restrict__ xtmp = (__v4df *__restrict)tmp;
    // hlinear_buffer<__v4df, 128, MP_IN> xxx(xtmp);
    // //dma_getn(tmp, ttmp, 512);
    // __v4df *__restrict__ ttmp = xxx.buffer;
    // __v4df *__restrict__ trsqrt = srsqrt;
    // memcpy(ttmp, tmp, 512*sizeof(__v4df));
    // for (int i = 0; i < 512; i += 16){
    //   trsqrt[i] = ttmp[i];
    //   trsqrt[i+1] = ttmp[i+1];
    //   trsqrt[i+2] = ttmp[i+2];
    //   trsqrt[i+3] = ttmp[i+3];
    //   trsqrt[i+4] = ttmp[i+4];
    //   trsqrt[i+5] = ttmp[i+5];
    //   trsqrt[i+6] = ttmp[i+6];
    //   trsqrt[i+7] = ttmp[i+7];
    // }
    // dma_putn(rsqrt, trsqrt, 512);
    // __v4df *__restrict__ ttmp = (__v4df *__restrict__)tmp;
    // __v4df *__restrict__ trsqrt = (__v4df *__restrict__)rsqrt;

    // hdistribute_iterate<128, 4>(0, 10000, hlinear_in(tmp), hlinear_out(rsqrt)).run<8>([&](int i) __cpe__ {
    //   *(__v4df*)(rsqrt + i) = *(__v4df*)(tmp + i);
    // });
    printf("1");
    _hiterate<1024, 4, 256*64>(0, n, hlinear_in(tmp), hlinear_out(rsqrt), [&](int i) __cpe__{
      *(__v4df*)(rsqrt + i) = v_invsqrt(*(__v4df*)(tmp + i));
    });
      //*(trsqrt + i) = v_invsqrt(*(ttmp + i));
  });

  double min = 0, max = 0;
  for (int i = 1; i < 100000; i ++){
    double check = tmp[i] * rsqrt[i] * rsqrt[i] - 1;
    if (check < min) {
    //  printf("%d %e\n", i, check);
      min = check;
    }
    if (check > max) {
    //  printf("%d %e\n", i, check);
      max = check;
    }
  }
 // printf("%e %e\n", min, max);
  return 0;
}
